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Abstract 


This is the abstract. 

Note: This is a working paper which will be expanded/updated frequently. All suggestions 
for improvement are welcome. The directory deleeuwpdx.net/pubfolders/gi6Evals has a pdf 
version, the bib 61c, the complete Rmd 61e with the code chunks, and the R source code. 

1 Introduction 

In a gihAnalysis (De Leeuw (2019)) we have n of observations on each of m variables. 
Each variable j has associated with it a (column-centered) basis matrix Hj and a closed 
polyhedral convex cone /C, in the space spanned by the columns of Hj. 

For each of the variables we want- to 6nd a vector of coefficients z 3 such that HjZj G /C y 
and such that the correlation matrix R of the transformed variables Uj = HjZj is, in 
some sense, optimal. 

To actually dehne optimality in the Gi6 context we associate with each gihAnalysis a 
dimension p and a partition of the set [M] = (1, 2, • • • , m} indexing the variables into r 
subsets Xi i, • • • , A4 r . In a traditional gihAnalysis we minimize gihLoss, dehned as 

a(X, Z,A)=± tr (A' - £ H j z, o')'(A - £ H^). (1) 

t =1 j&Mt j&Mt 

over X, over the cp, and over the Zj satisfying the cone restrictions (De Leeuw (2019)). The 
algorithm is a combination of alternating least squares and majorization (a.k.a MM). The 
Gih package in R implements this algorithm, for various choices of partitionings and bases. 

2 Gifi Eigenvalues 

In order to stay away from the trivial solution X = 0 and ZjCi'j = 0, which of course gives a 
zero loss, we need some form of normalization. 
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To discuss normalizations we use some additional notation. Define by horizontally 
concatenating the Hj with j G Ait- Define Z t as the direct sum of the z 3 with j G Ait. 
Thus Z t is block-diagonal, with the column-vectors Zj in the diagonal blocks. Also define 
Gt = ' HtZ t . Finally At vertically concatenates the row vectors % with j G Ai t . Thus loss 
function (1) becomes 

<?(X, Z,A) = ± tr (X - Q t A t )\X - g t At), (2) 

t =i 

De Leeuw (2019) discusses three different normalizations. The first one requires A^'A" = I but 
leaves A unnormalized, the second one requires YAt=\ A! t G' t GtAt = I and leaves X unnormalized. 
The third one requires both X and A to be normalized. All three normalizations give essentially 
the same solution for X, A, and Z. The three solutions only differ in the scaling of the p 
dimensions. 

2.1 Normalize X 

In the original Gift implementation, and in the R implementation of the homals package (De 
Leeuw and Mair (2009)), we normalize A" and leave A free. So we study this normalization 
first. 

The minimum of loss (2) over A is 

T 

miner (A", Z, A) = pT — ^ tr X'G t G+X, 

A t=l 

with Gt the Moore-Penrose inverse of Gt- Define Vt = G t Gti and define V+ as the average of 
the r matrices Vt- Then 

min miner(X, Z, A) = r(p - max tr X’V+X) = r(p - ^ A t (P*)), 

y(.'yv—1 A. y\.'y\.—1 , 

S= 1 

where the A S (P*) are the eigenvalues of V+, ordered from large to small. 

Of course we should not forget that the Gt, and thus T± and its eigenvalues, are all functions 
of the Zj. Thus the problem that still must be solved, after projecting out A and X from the 
loss function, is to maximize the sum of the p largest eigenvalues of V * over the Zj. 

2.1.1 Properties 

Since V t projects on the subspace C t spanned by the columns of Gt, it is idempotent and its 
eigenvalues are between zero and one. Thus the eigenvalues of P* are also between zero and 
one. An eigenvalue is equal to one if and only if its eigenvector is in the intersection of the 
Ct. In the terminology of Gifi (1990) that is if and only if it is in the meet of the Ct. Also 

0 < min miner (X.Z.A) < rp, 

where loss is equal to zero if and only if the C t have a p-meet, i.e. their intersection has 
dimension larger than or equal to p. 
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2.2 Normalize A 


One possible disadvantage of the original Gifi formulation is that the iterative minimization of 
loss has to be carried out in a space with dimension of order rip, the number of observations 
times the number of dimensions. This can, of course, be a very large number. Partly for that 
reason Tijssen and De Leeuw (1989) investigated minimization of Gifi loss using normalization 
of A, while leaving X free. 

First, observe that 

min a(X, Z,A) = J2 tr A t G' t GtAt - - tr A 8 Q' s Q t A t . 

t= 1 T t= 1 S=1 

Suppose C is the block matrix with blocks G' s Gt and T> is the block-diagonal matrix with 
blocks V t . If II is the m x m incidence matrix of the partition, i.e. n 3 c = 1 if two variables are 
in the same subset and xje = 0 otherwise, then T> is the elementwise or Hadamard product 
II x C. With these definitions 

min a(X, Z, A) = tr A!VA -tr A'CA, 

x v ' t 

and 

v 

^ nnn^niin cr(X, Z, A) = rp - J2 A(C, tV), 

S = 1 

A giftAnalysis find the Zj in their cones IC 3 such that the sum of the p largest 
eigenvalues of the generalized eigenvalue problem for the pair (.R , S), with S = ILvR, 
is maximized. 

2.2.1 Properties 

The eigenvalue problems defining a gifiAnalysis have some simple properties that are true for 
any partitioning of the variables. 

First, note that any gifiAnalysis is scale-free because for every diagonal D the eigenvalues 
of the pair ( DRD , DSD ) are the same as those of (R, S ). Thus, without loss of generality, 
we will restrict our covariance matrices to be correlation matrices. 

Second, the eigenvalues of the pair ( R , S) with S = II * R have simple upper bounds. 

Result 1: [Upper Bound]: If the partition II has r sets then 0 < R < rS. 

Proof: Partition the m columns of H using II, creating submatrices Hi, ■ • • , H r . Define the 
r linear combinations t q = H q y q and their average U. Then 

J2(tq - U )\t q - U) = y'Sy - r~ l y'Ry > 0. (3) 

q =1 
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Result 2: [Unity] The sum of the p largest eigenvalues of R and S = II ★ R is less than or 
equal to rp, with equality if and only if the r subspaces spanned by the transformed variables 
in the r sets have a /^-dimensional intersection (or meet). 

Proof: There is equality in (3) if and only if all t q are equal. Thus we have Ry = rSy if and 
only if the r subspaces have a non-trivial intersection. ■ 

2.2.2 Aspect Theory 

The eigenvalues of the pair (R,U-k R) are functions of the z r Using the terminology used in 
De Leeuw (1988c) the sum of the p largest eigenvalues is an aspect of the correlation matrix 
R. A gifiAnalysis chooses the Zj to maximize that particular aspect. 

A general majorization algorithm for aspects is outlined in De Leeuw (1988c), and implemented 
in R by Mair and De Leeuw (2010). It is proved to be convergent if the aspect is a convex 
functions of the correlation matrix. 

If the aspect / is convex then for any two correlation matrices R and R the majorization 
inequality is 

f(R) >/(i?) + tr (R-R)Vf(R). 

Suppose R is the current best estimate of the optimum correlation. If we find an R such that 
tr RVf(R) > tr RVf(R) then f(R) > f(R). What interests us in this paper is if the sum of 
the p largest eigenvalues of ( R , II * R) is a convex function of R, or, more precisely, for which 
partitions, if any, it is a convex function of R. 

It is well known that the sum of the p largest eigenvalues of a real symmetric matrix of 
order n is a convex function of its matrix argument. This follows easily from the Ky Fan 
representation of that sum as rnaxx tr X'RX, where X varies over the n x p orthonormal 
matrices. Indeed, tr X'RX is linear in R, and thus the sum of the largest eigenvalues is the 
maximum of a family of linear functions, which implies it is convex. The function we are 
optimizing is convex for the partition with exactly one variable in each set. We will study 
what happens for general partitions, for which the situation is less simple. 


3 Derivatives 

We can study convexity of Gifi’s eigenvalues by computing their first and second derivatives, 
assuming throughout that the eigenvalues we are differentiating are simple. Even if there is 
no convexity these derivatives are useful for optimization, not only of our sum-of-eigenvalues 
function, but for other functions of the eigenvalues used in multivariate analysis, such as the 
ones defined by Horst (1961) or Kettenring (1971). 

The correlation matrix R is a function of the coefficients z r We give the formulas for the 
derivatives of the aspect with respect to the coefficients in two steps, using the chain rule. 
First there is the derivative of the aspect as a function of the correlation matrix, and then 
there are the derivatives of the correlation matrix with respect to the coefficients. 
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We start with a slightly more general formulation than is needed for a gihAnalysis. Consider 
the problem Ay = A By and y' By = 1 where 

I< 

A(9) = A 0 + E 9kAk, 

k=i 

K 

B{9) = B 0 + E QhBk- 

k =1 

Here Ap. and Bp. are symmetric matrices. Thus both A and B are defined as functions of the 
vector with the K elements 6p. This is more general than what is needed for Gih, because in a 
gihAnalysis A and B have very specific forms. We will specialize our results after deriving the 
more general ones, which are useful in other contexts (for example when deriving confidence 
regions for correspondence analysis solutions, see De Leeuw (2007)). 

For the first and second derivatives we use the formulas in De Leeuw (2007), although the 
results given there also appear in hundreds of other places in the literature of functional 
analysis, numerical analysis, multivariate statistical analysis, optimization, and engineering. 
The approach in De Leeuw (2007) is similar, for example, to that of Friswell (1994), who also 
discusses how to generalize the formulas to higher order derivatives. Shapiro (1985) shows 
how to derive the required results directly from the implicit function theorem. 

As an aside, there are also numerous results available, both in the mathematical and 
the engineering literature, for the case of multiple eigenvalues, which are generally only 
directionally differentiable functions of their matrix argument. For eigenvalues of symmetric 
matrices formulas for the second-order directional derivatives have been given by Torki (2001) 
and Zhang, Zhang, and Xiao (2013). We will not use these results here, and just keep them 
in mind for possibly later extensions. 

First some definitions. Define the K vectors (Ap — X s (9)Bp)y s (9) of length n, and concatenate 
them vertically to the K x n matrix H s (9). Also define the /7-element vector g s (9) with 
elements y' s (9)Bpy s (9) and the //-element vector f s (9) = H s (9)y s {9). Finally, again following 


De Leeuw (2007), we define a type of generalized inverse 

W s {9) = ( A(9) - A s {9)B{9))~ = Y(9)(K(9) - A S (9)I) + Y'(9), (4) 

where Y(9) and A (9) are complete sets of eigenvectors and eigenvalues for the pair (A(9), B(9)), 
and superscript + is the Moore-Penrose inverse. 

With these definitions we have nice compact matrix expressions for the derivatives. 

V\ s (9) = H s (9)y s (9) = f s (9), (5) 

Vy s (9) = -W S {9)H S {6) - \y s {0)g' s {6 ), (6) 

V 2 X s (9) = -2 H S {9)W S {9)H' S {9) - g s {9)f' s {9) - f s {9)g' s {9). (7) 

3.1 Signature 


The signature (sometimes called the inertia) of a real symmetric matrix is the triple (n, u, 6) 
with the number of positive, negative, and zero eigenvalues. Sylvester’s Law of Inertia says 
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that for any square non-singular X the signature of X'AX is equal to the signature of A. In 
Dancis (1986) this is generalized to arbitrary X, not necessarily square or non-singular. We 
give a different proof of their theorem 3.1. 

Result 3: [Poincare] Suppose X is an n x m matrix of rank r and A is a symmetric 
matrix of order n. If the signature of A is (n, v, S) and the signature of X'AX is ft, u, 6 then 
7T — (n — r) < 7T < 7T and v — (n — r) < v < v. 

Proof: Suppose the singular value decomposition of X is 


X 


K 


A 0 
0 0 


L', 


with K orthonormal of order n, L orthonormal of order m, and A diagonal and positive 
definite of order r. Then 


X'AX 


AK[AK\A 0 
0 0 


L\ 


where A\ is n x r with K\ K x = I. By applying Sylvester’s Law of Inertia the signature 
of X'AX is the signature of AK[AK 1 A plus (0,0, m — r). Applying Sylvester’s Law again 
the signature of AK[AK 1 A is the signature of K[AK 1 . Now apply Poincare’s Separation 
Theorem (see, for example, Abadir and Magnus (2005), exercise 12.46) to deduce that the 
eigenvalues AT > • • • > fi r of K[AK 1 are related to the eigenvalues Ai > • • • > \ n of A by the 
inequalities A n -r+i < AA A*. Since Hi A A, it is clear that n < n. Thus at > 0 if A n -r+i > 0 
if n — r + i < n or % < n — (n — r). Thus n > n — (n — r). The result for v follows by applying 
the same reasoning to —A, which of course has signature (u, n,5). ■ 


Result 4: Signature Suppose the K x (n + 1) matrix | H s {9) ,g iS (0)j has rank n + 1. Then 
V 2 \ s {6) has n — s + 1 positive and s negative eigenvalues. 


Proof: We apply result 5 to the second derivatives V 2 X s (9). Write 


V 2 \ s (0) 


H a (6) f a (6) g s (6) 


'2 W a {6) 
0 
0 


0 0 
0 1 
1 0 


H' s (6) f'(9) g'M 


Now the signature of W s {6 ) is (s — 1, n — s, 1) and thus the signature of 

'2 W 8 {0) 0 O' 

0 0 1 
0 1 0 


is (s, n — s + 1,1). Taking the negative gives (n — s + 1, s, 1). The matrix 

\H S (0) f a (9) g,(0) 

is K x (n + 2), but since f s (9) is in the column space of H s {9 ) its rank is less than or 
equal to n + 1. We assume it is actually n + 1. Then result 5 gives us7r — l<7r<7ror 
n — s < n < n — s + 1. And v — l<z/<z/ors — 1 < i> < s. 
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Any vector z satisfying both z'H s (9) = 0 and z'g s (9 ) = 0 satisfies V 2 \ s {9)z = 0. Since we 
assume that r = n + 1 there are n + 1 non-zero eigenvalues It follows that tx = n — s + 1 and 
v — s. ■ 

It is abundantly clear from the results on signature that we can forget about- proving for 
egeneral partitions the convexity of the eigenvalues and their sums. 

3.2 Computation 

The R function computeABQ in the appendix computes and returns the first derivative of a 
single eigenvalue and corresponding eigenvector with respect to 9, using formulas (5), and (6). 
It also computes the second derivatives of the eigenvalue, using (7). In order to make sure 
that everything is OK, we also use checkABQ, which uses the numDeriv package (Gilbert and 
Varadhan (2019)) to compute and return the same quantities using numerical differentiation. 

We give the eigenvalues of the matrix of second derivatives for the partition (1,1, 2, 3,3, 3) 
and for the first three eigenvalues, using both computeABQ and checkABQ. 

rab <- makeAB (c(l,l,2,3,3,3)) 
h <- computeAB (rab$r, rab$a, rab$b, 1) 
mprint (eigen (h$hlf)lvalues, d = 6, w = 8, f = 


## 

[1] 

+22.434684 +17.508448 

+10.774439 

+7.727820 

+0.580109 

+0.078904 

## 

[7] 

+0.000000 +0.000000 

+0.000000 

+0.000000 

+0.000000 

- 0.000000 

## 

[13] 

- 0.000000 - 0.000000 

-0.511440 




h < 

- checkAB (rab|r, rab|a, rab|b, 1) 




mprint 

(eigen (h|hln)lvalues , 

d = 6, w = 

8, f = "+- 

") 


## 

[1] 

+22.434684 +17.508448 

+10.774439 

+7.727820 

+0.580109 

+0.078904 

## 

[7] 

+0.000000 +0.000000 

+0.000000 

- 0.000000 

- 0.000000 

- 0.000000 

## 

[13] 

- 0.000000 - 0.000000 

-0.511440 





h <- computeAB (rab$r, rab$a, rab$b, 2) 
mprint (eigen (h$hlf)lvalues, d = 6, w = 8, f = 


## 

[1] 

+66.324525 +25.916113 

+7.638237 

+0.745604 

+0.572377 

+0.000000 

## 

[7] 

- 0.000000 - 0.000000 

- 0.000000 

- 0.000000 

- 0.000000 

- 0.000000 

## 

[13] 

-0.000000 -0.606543 

-18.251480 




h < 

- checkAB (rab|r, rab|a, rab|b, 2) 




mprint 

(eigen (h|hln)lvalues , 

d = 6, w = 

8, f = "+- 

") 


## 

[1] 

+66.324525 +25.916113 

+7.638237 

+0.745604 

+0.572377 

+0.000000 

## 

[7] 

+0.000000 +0.000000 

+0.000000 

+0.000000 

- 0.000000 

- 0.000000 

## 

[13] 

-0.000000 -0.606543 

-18.251479 





h <- computeAB (rab|r, rab|a, rab|b, 3) 
mprint (eigen (h|hlf)lvalues, d = 6, w = 8, f = 


7 


## [1] +59.742434 +29.991981 +10.385545 +0.043659 +0.000000 +0.000000 
## [7] +0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 
## [13] -0.044127 -19.285103 -62.029367 

h <- checkAB (rab$r, rab$a, rab$b, 3) 

mprint (eigen(h$hln)$values, d=6, w=8, f= "+-") 


## [1] +59.742434 +29.991981 +10.385545 +0.043658 +0.000000 +0.000000 
## [7] +0.000000 +0.000000 +0.000000 -0.000000 -0.000000 -0.000001 
## [13] -0.044128 -19.285103 -62.029367 

Clearly the results of computeABQ and checkAB() are basically the same, except for some 
minor numerical differences. Moreover the signatures agree with our theoretical results. 
For this small example the results using numerical derivatives are quite impressive. The 
computing the exact (formula) results is slightly more efficient, and we expect that using 
formulas will have more of an advantage in efficiency and precision for larger examples. 


3.3 Partials 

We now specialize our results to compute partial derivatives with respect to the lower-diagonal 
elements of the correlation matrix R = {pij}, with i > j. Thus we have K = \m(m — f), 
where m is the number of variables. We use double-subscripting for the K matrices A t]: 
which are m x m. They are of the form A tJ = e^e)- + e^-e), with e* and e 3 unit vectors. The Aij 
with i > j have elements (i, j) and (j,i) equal to +1, while all other elements are zero. Thus 

A(p) — I + Pi j A^, 

l<jf<2<m 

5(p) I ~ j- ^ ^ • 


There is no need to compute and store A^ and D l3 in this case, we can simply substitute their 
functional form into equations (5), (6) and (7). This obviously saves both multiplications and 
memory. 

Thus, for i > j, the derivatives are 


T^(ij)A s (p) — 2(1 \ s (p)i7ij)y s i(p)y s j(p), (8) 

and 

v {H)y(.s{p) = -(1 - X s (p)TTij)(yjs{p)w si t(p) + y is (p)w sje (p)) - -K ij y is (p)y is (p)y js (p), (9) 

where W s (p) is defined by (4). 

The second derivatives are 

v mk iMp) — 

‘2'D(kl)X s (p)'H'ijyis(p)yjs(p) T 2(1 A s (p)7Tjj ) (©(fcql/js (p)?/j S (p) + yis(p)'D(kl)yjs(p)) (10) 


3.4 Computation 

In the appendix we give the R function computePQ which implements formulas (??), (??) 
and (??). We use a small example with four variables, in two sets, to compute first and 
second partials with both computeABQ and computePQ. 

p <- c(l, 1, 2, 2) 
rab <- makeAB(p) 

hab <- computeAB (rab|r, rab|a, rab|b, 1) 
hpa <- computeP (rab|r, p, 1) 

We now compare the two results to check if our formulas and our computations are correct. 
The maximum absolute difference between the two results for the gradient of the eigenvalue is 
1.0408340856 x 10 -17 . For the gradient of the eigenvector it is 2.4980018054 x 10 -15 , and for 
the hessian it is 4.0523140399 x 10 -15 . Again, basically the same, so things seem to be OK. 

3.5 Special Cases 

3.5.1 Component Analysis 

If 7 Tij = 0 for all % ^ j (as in the case of nonlinear principal component analysis) 

V 2 \ s (6) = -2H'(0)(A(0) - A s (0)I) + H s (0) 

It follows, assuming that H s {9 ) is of rank n, that the signature of T> 2 \ s {6) is (n — s,s — 
1, |[n — 1 )(n — 2)). Thus the largest eigenvalue is a convex function and the smallest is a 
concave function of the matrix argument. 

We verify the signature using computePQ on a Spearman correlation matrix of order six, for 
the first three eigenvalues. 

p <- 1:6 

z <- ( 1 : 6 ) / 10 

r <- outer (z, z) [outer(l: 6, 1:6, "<")] 
hi <- computeP (r, p, 1) 

mprint (eigen(hl$hlf) lvalues , d = 6, w = 8, f = 

## [1] +2.787716 +2.321184 +1.756067 +1.252285 +0.880824 +0.000000 +0.000000 

## [ 8 ] + 0.000000 + 0.000000 + 0.000000 - 0.000000 - 0.000000 - 0.000000 - 0.000000 
## [15] -0.000000 

h2 <- computeP (r, p, 2) 

mprint (eigen(h2$hlf) lvalues , d = 6, w = 8, f = "+-") 

## [1] +51.743207 +20.223939 +10.841689 +6.722973 +0.000000 +0.000000 

## [7] +0.000000 +0.000000 +0.000000 -0.000000 -0.000000 -0.000000 

## [13] -0.000000 -0.000000 -2.782146 
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h3 <- computeP (r, p, 3) 

mprint (eigen(h3$hlf) lvalues , d = 6, w = 8, f = 


## 

[1] 

+27.249588 

+12.931445 

+7.521237 

+0.000000 

+0.000000 

+0.000000 

## 

[7] 

+0.000000 

+0.000000 

+0.000000 

-0.000000 

-0.000000 

-0.000000 

## 

[13] 

-0.000000 

-2.247776 

-51.689503 




We 

can also check that the sum of the first two and the sum 

of the first three eigenvalues is 

convex. 







mprint 

(eigen(hl$hlf+h2$hlf)$values, d = 

6, w = 8, 

f = "+-») 


## 

[1] 

+51.763142 

+20.233969 

+10.848115 

+6.727161 

+2.305156 

+1.747129 

## 

[7] 

+1.246273 

+0.876793 

+0.000000 

+0.000000 

+0.000000 

-0.000000 

## 

[13] 

-0.000000 

-0.000000 

-0.000000 




mprint 

(eigen(hl$hlf+h2$hlf+h3$hlf)$values, d = 6, 

w = 8, f = 

"+-") 

## 

[1] 

+27.518920 

+20.165318 

+12.995857 

+10.841753 +7.548198 

+6.725746 

## 

[7] 

+1.648177 

+1.211152 

+0.857608 

+0.000000 

+0.000000 

+0.000000 

## 

[13] 

-0.000000 

-0.000000 

-0.000000 





3.5.2 Regression 

Suppose there are two sets and the first set contains only a single variable (as in the case of 
nonlinear multiple regression). We can partition the correlation matrix as 

1 r' 
r R 


The largest eigenvalue is 1 + a /r'Rr, the smallest 1 — Vr'RrW, and there are n — 2 eigenvalues 
equal to 1. Now 1 + \f r'Rr is not a convex function of the correlation matrix, but r'R~ x r is 
convex (De Leeuw (1988c)). It follows that (1 — Ai(6*)) 2 is convex. Obviously maximizing the 
largest eigenvalue is the same as maximizing the aspect r'R~ 1 r. 


Our usual result on signatures applies in the regression case as well. 


p <- 1:6 

z <- (1 : 6) / 10 

r <- outer (z, z) + diag(l - z ~ 2) 
mprint (r) 


## 

## [1J 
## [ 2 ,] 
## [3,] 
## [4,] 
## [5,] 
## [ 6 ,] 


[, 1 ] 

+ 1.000000 

+ 0.020000 

+0.030000 

+0.040000 

+0.050000 

+0.060000 


[, 2 ] 

+ 0.020000 

+ 1.000000 

+0.060000 

+0.080000 

+ 0.100000 

+ 0.120000 


[ ,3] 

+0.030000 

+0.060000 

+ 1.000000 

+ 0.120000 

+0.150000 

+0.180000 


[>4] 

+0.040000 

+0.080000 

+ 0.120000 

+ 1.000000 

+ 0.200000 

+0.240000 


[, 5] 

+0.050000 

+ 0.100000 

+0.150000 

+ 0.200000 

+ 1.000000 

+0.300000 


[, 6 ] 

+0.060000 

+ 0.120000 

+0.180000 

+0.240000 

+0.300000 

+ 1.000000 
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lr <- eigen (r)lvalues[1] 
r <- r[outer(l:6, 1:6, "<")] 
p <- c(l,2,2,2,2,2) 
hp <- computeP (r, p, 1) 

mprint (eigen(hp$hlf)lvalues, d = 6, w = 8, f = 

## [1] +19.356464 +16.685695 +15.117584 +14.179460 +0.290875 +0.000000 

## [7] +0.000000 +0.000000 +0.000000 +0.000000 +0.000000 +0.000000 

## [13] -0.000000 -0.000000 -0.262528 

hq <- computeP (r, p, 2) 
mprint (hqlglf) 

## [ 1 ] - 0.000000 - 0.000000 - 0.000000 + 0.000000 + 0.000000 + 0.000000 - 0.000000 
## [ 8 ] - 0.000000 - 0.000000 + 0.000000 + 0.000000 + 0.000000 + 0.000000 - 0.000000 
## [15] +0.000000 

We expect convexity for p{9) = (1 — A^#)) 2 . Now 

v 2 p(d) = -2(1 - A 1 (0))D 2 A 1 (0) + 2V\ 1 (9)(V\ 1 (9))', 


and we can thus use the results of computeAB() to verify convexity in the regression case. 

hlmp <- 2 * ( 1 - lr) * hplhlf + 2 * outer (hplglf, hplglf) 

mprint (eigen (hlmp)lvalues, d = 6, w = 8, f = "+-") 

## [1] +26.069560 +22.469558 +20.355243 +19.090233 +1.353669 +0.000000 
## [7] +0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 
## [13] -0.000000 -0.000000 -0.090875 

mab <- makeAB(p) 

hab <- computeAB(r, mab|a, mab|b, 1) 
print(max(abs (hplglf - hablglf))) 

## [1] 6.938893904e-18 

print(max(abs (hplgyf - hablgyf))) 

## [1] 3.330669074e-16 

print(max(abs (hplhlf - hablhlf))) 

## [1] 4.9960036lie-16 

hlmab <- -2 * ( 1 - lr) * hablhlf + 2 * outer (hablglf, hablglf) 
print(max(abs (hlmp - hlmab))) 

## [1] 7.21644966e-16 
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mprint (eigen (hlmab) lvalues , d = 6, w = 8, f = 

## [1] +26.069560 +22.469558 +20.355243 +19.090233 +1.353669 +0.000000 
## [7] +0.000000 +0.000000 +0.000000 +0.000000 -0.000000 -0.000000 
## [13] -0.000000 -0.000000 -0.090875 


4 GifiA Derivatives 


In a gifiAnalysis 


Pij(z) = 


zfiijZj 




where the C, :! = //'// y are submatrices of the Burt matrix (remember the Hi are column- 
centered). By the chain rule 


V k \ s (z) = T)(i j )X s (p(z))V k p ij (z). 

l<j‘<i<n 


The function computeZ() combines (11) with 

Pi j (^) 


y/ z i D i Zi 


{I ~ DiZityCijZj, 


v jpij( z ) = r— ( J - D^fyCjiZi ,. 
\J Z 3 D 3 Z 3 


( 11 ) 


( 12 ) 

(13) 


This is not the most efficient way of computing the partials, but it gets the job done. In the 
same way, by the chain rule, 


H k y s (z) = J2 V {ij) y s {p(z))V k pij(z). (14) 

l<j‘<i<n 


For the second partials we have, again using the chain rule, 


^fc£-^s(^) } ) ^ ) ^(n)(vo)^s{p{z'))H k pij[z)'D uPpq^z) -\- ^ A (p(z)\H k ipjj (/^). (15) 

j<iq<p j<i 


Once again, the computeZ() function does not aim to be efficient. It directly implements the 
chain rule $(15), using the partials from computePQ, and in addition 


'DijPij (*) 


'^iiPij (4 


DjjPij ( z ) 


4 A z i 




{I - D i z i z 'i)Cij(I ~ ZjZjDj), 


z'i.Di Zi 

1 

Z 3 D 3 Z 3 


3 pijHiZi^iDi ( D,jz^ZjCjj T CjjZjzjDj ) pijDi , 

3 pijDjZjZjDj — (DjZjZiCij + CjiZiZjDj) — PijDj 


(16) 

(17) 

(18) 
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4.1 Bilinearizability 

It is interesting that if there are Zj with z'jDjZj = 1 that satisfy the simultaneous bilinear 
regression conditions CkjZj = r k jD k z k (see De Leeuw (1988b) and De Leeuw (1988a)) then by 
(9) and (10) V k pij{z ) = 0 and thus V k \ s (z) = 0 and T> k y s (z), no matter what the partition 
is. It follows that a system of bilinearizing Zj corresponds with a stationary solution for any 
aspect that is a function of the correlation coefficients. 

Not necessarily a local maximum of minimum, of course. In fact from (16), (17) and (18) we 
see that 


Djj f)jj (z'j Cij f)jj (z)Dj Z-I Zj Dj , 

^iiPij (^) Pij (t) D{ZiZj^D j), 

Pij (t) Pij ( z ) ( Dj Dj Zj Zj Dj ). 

Thus, for simultaneously bilinearizing scores z v 

D ke X s (z) = 2 5Z(1 - K{p{z))TT i j)y si (p(z))y sj (p(z))T> ke p ij (z). (19) 

3<i 


4.2 Computation 

icat <- c (2, 2, 4, 4) 

p <- c (1, 1, 2, 2) 

hz <- makeHZ (100, 4, icat) 

cf <- computeZ (hz$z, hz$h, hz$r, p, 1) 

cn <- checkZ (hz$z, hz$h, p, 1) 

print(max(abs(uniist (cf$glf) - cn$gln))) 

## [1] 1.473018209e-10 

print(max(abs(matList(cf$gyf ) - cn$gyn))) 

## [1] 0.8471463466 

hi <- ifelse (outer(sample (1:4, 100, replace = TRUE), 1:4, "=="), 1, 0) 

h2 <- ifelse (outer(sample (1:6, 100, replace = TRUE), 1:6, "=="), 1, 0) 

hi <- apply (hi, 2, function (x) x - mean (x)) 

h2 <- apply (h2, 2, function (x) x - mean (x)) 

hi <- svd (hi) $u [, -4, drop = FALSE] 

h2 <- svd (h2)$u[, -6, drop = FALSE] 

cl2 <- crossprod (hi, h2) 

s <- svd (cl2) 

zl <- s$u[,l] 

z2 <- s$v[,l] 

h <- list (hi, h2) 

z <- list (zl, z2) 
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r <- zl 1*1 cl2 1*1 z2 
mprint (r) 

## [, 1 ] 

## [1,] +0.238971 

p <- c (1 , 2) 

zcompu <- computeZ(z, h, r, p, 1) 
hess <- matrix(0, 8, 8) 
hess [1:3,1: 3] <-zcompu$hlf[ [1] ] 
hess [4 : 8 ,1: 3] <-zcompu$hlf[ [2] ] 
hess [4 : 8 , 4: 8] <-zcompu$hlf[ [3] ] 
hess[1 : 3 ,4: 8] <-t (zcompu$hlf[ [2] ]) 
mprint (eigen(hess) $values) 

## [1] -0.000000 -0.000000 -0.043703 -0.139061 -0.238971 -0.238971 -0.338882 
## [8] -0.434239 

zcompu <- computeZ(z, h, r, p, 2) 
hess [1:3,1:3] <-zcompu$hlf[[1]] 
hess [4 : 8 ,1: 3] <-zcompu$hlf[ [2] ] 
hess [4 : 8 , 4: 8] <-zcompu$hlf[ [3] ] 
hess[1 : 3 ,4:8] <-t (zcompu$hlf[ [2] ]) 
mprint(eigen (hess)$values) 

## [1] +0.434239 +0.338882 +0.238971 +0.238971 +0.139061 +0.043703 +0.000000 
## [ 8 ] + 0.000000 


5 Appendix: Duality Gap 


The result in this appendix may not be of much practical use, but it’s a nice example of the 
duality gap, so I’ll include it anyway. The primal problem is 

max tr X'AX = max min tr X'AX — tr M(X'BX — I). 

X'BX=I X M£S n v 7 

We now study the Lagrange dual problem 

min maxtr X'AX — tr M(X'BX — I). 

M£S n X K 1 

Using Kronecker products this becomes 


maxtr x'[(I (g) A) — (M (g) B)\x + tr M 


tr M if M <g) B > I <g) A, 
Too otherwise. 


Now 


min 

MeS n 


{tr M | M <g B > I ® A} = pXpiB^A), 
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and thus 


min max tr X'AX - tr MiX'BX - I) < max min tr X'AX - tr MiX'BX - I), 

Mes n x K ' x MeS n 

with equality if and only if the p largest eigenvalues of B~ Y A are all equal. 


6 Appendix: MaxMax Duality 


There is another result which I am just reporting here because I have it sitting around. First 


tr X'AX = max 2 tr Y'AX - tr Y'AY. 

Y 


Thus 

max tr X'AX = max max 2 tr Y'AX - tr Y'AY = max max 2 tr Y'AX - tr Y'AY 

X'BX=I X'BX=I Y Y X'BX=I 

Now 

max 2 tr Y’AX - tr Y'AY = 2 tr (Y'AB~ l AY)^ - tr Y'AY. 

and thus 

max tr X'AX = max 2 tr (Y'AB^AY)^ - tr Y’AY. 

X’BX=I Y 

Note that tr (Y'AB~ 1 AY)^ is the sum of the singular values, i.e. the trace norm, of Y’AB~^. 

Thus we have replaced a constrained problem with an, admittedly more complex, uncon¬ 
strained problem. 

7 Appendix: Code 

7.1 auxiliary. R 

library (numDeriv) 

repList <- function(x, n) { 
z <- listO 
for (i in l:n) 

z <- c(z, list(x)) 
return(z) 

} 

sumList <- function (x, a) { 
m <- length (x) 
nr <- nrow(a[[l]]) 
nc <- ncol(a[ [1] ]) 
sa <- matrix (0, nr, nc) 
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for (j in l:m) { 

sa <- sa + x[j] * a[[j]] 

} 

return (sa) 


matList <- function (x) { 
m <- length (x) 
n <- nrow (x [ [1]] ) 
a <- matrix(0 ,n, 0) 
for (j in l:m) { 

a <- cbind (a, x[ [j]]) 

} 

return (a) 

} 

psplit <- function (z, s) { 

g <- 1 + rowSums (outer (1 : length(z) , cumsum(s) , ">")) 
h <- split (z, g) 
names (h) <- NULL 
return (h) 


geig <- function (a, b) { 
s <- solve (chol (b)) 
h <- crossprod (s, a %*% s) 
e <- eigen (h) 
y <- s e$vectors 

return (list (values = e$values, vectors = y)) 

} 

mprint <- function (x, d = 6, w = 8, f= { 

print (noquote (formate ( 

x, di = d, wi = w, fo = "f", flag = f 
))) 


7.2 computeAB.R 

makeAB <- function (p) { 
m <- length (p) 

rmat <- cor (matrix (rnorm (1000 * m), 1000, m)) 
r <- rmat [outer (1:m, l:m, "<")] 


16 


a <- repList(matrix(0 , m, m), m * (m - 1) / 2) 
b <- repList(matrix (0, m, m), m * (m - 1) / 2) 
k <- 1 

for (i in 2:m) { 

for (j in 1: (i - 1)) { 
a[ [k] ] [i, j] <- 1 
a[ [k] ] [ j , i] <- 1 
if (p[i] == p [j]) { 
b [ [k] ] [i, j] <- 1 
b [ [k] ] [j , i] <- 1 

} 

k <- k + 1 

> 

} 

return (list (r = r, a = a, b = b)) 

} 

computeAB <- function (x, a, b, s) { 
m <- length (x) 
n <- nrow (a[[1]]) 
ax <- sumList (x, a) + diag (n) 

bx <- sumList (x, b) + diag (n) 

e <- geig (ax, bx) 

1 <- e$values 
y <- e$vectors 
ys <- y[, s, drop = TRUE] 
v <- ifelse (s == l:n, 0, 1 / (1 - l[s])) 

w <- tcrossprod (y 1*1 diag (v), y) 

hml <- sapply(a, function (x) 
x 1*1 ys) 

hm2 <- sapply(b, function (x) 
x 1*1 ys) 

hmat <- hml - 1[s] * hm2 
gvec <- sapply (b, function (x) 
drop (ys 1*1 x 1*1 ys)) 
gyf <- -w 1*1 hmat - outer (ys, gvec) / 2 
gif <- drop (ys 1*1 hmat) 
hlf <- 

-2 * crossprod (hmat, w 1*1 hmat) - outer (gif, gvec) - outer (gvec, gif) 
return (list ( 
gif = gif, 
gyf = gyf, 
hlf = hlf 

)) 
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> 


checkAB <- function (x, a, b, s) { 
fl <- function (x, a, b, s) { 
n <- nrow(a[[l]]) 
ax <- sumList (x, a) + diag (n) 

bx <- sumList (x, b) + diag (n) 

e <- geig (ax, bx) 
return (e$values[s]) 

} 

fy <- function (x, a, b, s) { 
n <- nrow(a[[l]]) 
ax <- sumList (x, a) + diag (n) 

bx <- sumList (x, b) + diag (n) 

e <- geig (ax, bx) 

k <- which.max (abs (e$vectors[, s])) 

return (sign (e$vectors[k, s]) * e$vectors[, s]) 

} 

gin <- grad ( 

fl, 

x, 

side = NULL, 
a = a, 
b = b, 
s = s 

) 

gyn <- jacobian ( 

fy, 

x, 

side = NULL, 
a = a, 
b = b, 
s = s 

) 

hln <- hessian (fl, x, a = a, b = b, s = s) 
return (list (gin = gin, 
gyn = gyn, 
hln = hln)) 


7.3 computeP.R 
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computeP <- function (r, p, s) { 
m <- length (p) 
mm <- m * (m-1) / 2 
a <- matrix (0, m, m) 
a[outer(l:m, l:m, "<")] <- r 
a <- a + t(a) 
diag(a) <- 1 

pp <- ifelse(outer(p, p, "=="), 1, 0) 
b <- pp * a 
e <- geig (a, b) 

1 <- e$values 

y <- e$vectors 

ys <- y[, s, drop = TRUE] 

Is <- 1 [s] 

v <- ifelse (s == l:m, 0, 1 / (1 - Is)) 

w <- tcrossprod (y °/ 0 *7o diag (v), y) 

gg <- rep (0, mm) 

gif <- rep (0, mm) 

gyf <- matrix (0, m, mm) 

hlf <- matrix (0, mm, mm) 

k <- 1 

for (i in 2:m) { 

for (j in 1: (i - 1)) { 

gif [k] <- 2 * (1 - pp[i, j] * Is) * ys[i] * ys [j] 
aux <- 

-(1 - pp [i, j] * Is) * ((ys [j] * w[, i]) + (ys [i] * w[, j])) 
gyf [, k] <- aux - pp[i, j] * ys [i] * ys[j] * ys 
1 <- 1 

for (u in 2:m) { 

for (v in 1: (u - 1) ) { 

} 

} 

k <- k + 1 

> 

} 

k <- 1 

for (i in 2:m) { 

for (j in 1: (i - 1)) { 

facl <- -2 * pp[i, j] * ys [i] * ys[j] 
fac2 <- 2 * (1 - pp[i, j] * Is) 

1 <- 1 

for (u in 2:m) { 

for (v in 1: (u - 1) ) { 
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fac3 <- (ys [j] * gyf[i, 1]) + (ys [i] * gyf 
hlf[k, 1] <- facl * gif [1] + fac2 * fac3 
1 <- 1 + 1 

} 

} 

k <- k + 1 

> 

} 

return (list (gif = gif, 
gyf = gyf, 
hlf = hlf)) 


checkP <- function (r, p, s) { 
fl <- function (r, p, s) { 
m <- length (p) 
mm <- m * (m - 1) / 2 
a <- matrix (0, m, m) 
a[outer(l:m, l:m, "<")] <- r 
a <- a + t(a) 
diag(a) <- 1 

pp <- ifelse(outer (p, p, "=="), 1, 0) 
b <- pp * a 
e <- geig (a, b) 
return (e$values[s]) 

} 

fy <- function (r, p, s) { 
m <- length (p) 
mm <- m * (m - 1) / 2 
a <- matrix (0, m, m) 
a[outer(l:m, l:m, "<")] <- r 
a <- a + t(a) 
diag(a) <- 1 

pp <- ifelse(outer (p, p, "=="), 1, 0) 
b <- pp * a 
e <- geig (a, b) 

k <- which.max (abs (e$vectors[, s])) 

return (sign (e$vectors[k, s]) * e$vectors[, s]) 

} 

gin <- grad ( 

fl, 

x, 

side = NULL, 
a = a, 


j, 1]) 
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b = b, 
s = s 

) 

gyn <- jacobian ( 

fy, 

X, 

side = NULL, 
a = a, 
b = b, 
s = s 

) 

hln <- hessian (fl, x, a = a, b = b, s = s) 
return (list (gin = gin, 
gyn = gyn, 
hln = hln)) 


7.4 computeZ.R 

makeHZ <- function (n, m, k) { 
h <- as.list (1:m) 
z <- as.list (l:m) 
if (length(k) == 1) { 
k <- rep (k, m) 

} 

hzz <- matrix (0, n, m) 
for (j in l:m) { 
hh <- 

if else (outer (sample ( 1: k[j], n, replace = TRUE), l:k[j], "=="), 1, 0) 
mh <- apply (hh, 2, mean) 
hh <- hh - outer (rep(l, n), mh) 
hh <- svd (hh)$u[, — k [j ] , drop = FALSE] 
zz <- rnorm (k[j] - 1) 
hz <- drop (hh 0 /o*°/ 0 zz) 
zz <- zz / sqrt (sum (hz “2)) 
hzz[, j] <- hz 
h[[j]] <- hh 
z[[j]] <- zz 

} 

r <- cor (hzz) [outer (1:m, l:m, "<")] 
return (list(h = h, z=z, r = r)) 

} 

ff <- function (z, h, p) { 
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m <- length (h) 
n <- nrow (h[ [1] ]) 
hz <- matrix (0, n, m) 
for (j in l:m) { 

hz[, j] <- h [ [j] ] 7. *7. z [ [j] ] / sqrt (sum (z [ [j] ] ~ 2)) 

} 

ax <- crossprod (hz) 

pp <- ifelse (outer (p, p, "=="), 1, 0) 
bx <- pp * ax 
e <- geig (ax, bx) 
return (e) 


computeZ <- function (z, h, r, p, s) { 
m <- length (z) 
cp <- computeP (r, p, s) 
hlf <- as.list (l:(m * (m + 1) / 2)) 
gif <- as.list (l:m) 
gyf <- as.list (l:m) 
kl <- 1 

for (k in l:m) { 

drk <- drmat (z, h, k) 

gif [ [k] ] <- drop (cp$glf */,*•/, drk) 

gyf[[k]] <- cp$gyf %*% drk 

for (1 in 1: k) { 

drl <- drmat (z, h, 1) 

hlf[[kl]] <- crossprod (drk, cp$hlf 7o*7* drl) 
dkl <- ddmat (z, h, k, 1) 

hlf[[kl]] <- hlf[[kl]] + sumList (cp$glf, dkl) 
kl <- kl + 1 

> 

} 

return (list (gif = gif, gyf = gyf, hlf = hlf)) 

> 

drmat <- function (z, h, k) { 
m <- length (z) 
mm<-m* (m-1) / 2 
dr <- matrix (0, mm, length (z[[k]])) 

1 <- 1 

for (i in 2:m) { 

for (j in 1: (i - 1)) { 
if (i == k) { 

dr[1,] <- dirij (z, h, i, j) 


22 


} 

if (j == k) { 

dr[1,] <- djrij (z, h, i, j) 

} 

1 <- 1 + 1 


} 

return (dr) 


ddmat <- function (z, h, k, 1) { 
m <- length (z) 
inm<-m* (m-1) / 2 

dd <- repList(matrix (0, length (z[[k]]), length(z [[1]])), mm) 
kl <- 1 

for (i in 2:m) { 

for (j in 1: (i - 1)) { 

if ((i == k) && (j == 1)) { 

dd[[kl]] <- dijrij (z, h, i, j) 

} 

if ((i == k) && (i == 1)) { 

dd[[kl]] <- diirij (z, h, i, j) 

} 

if ((j == k) && (j == l)) { 

dd[[kl]] <- djjrij (z, h, i, j) 

} 

kl <- kl + 1 

> 

} 

return (dd) 

} 


dirij <- function (z, h, i, j) { 
cij <- crossprod (h[[i]], h[[j]]) 
cii <- crossprod (h[[i]]) 
czj <- drop (cij 1*1 z[[j]]) 

return (czj - sum (czj * z[[i]]) * drop (cii 1*1 z[[i]])) 

} 


djrij <- function (z, h, i, j) { 
cji <- crossprod (h[[j]], h[[i]]) 
cjj <- crossprod (h[[j]]) 
czi <- drop (cji 1*1 z[[i]]) 
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return (czi - sum (czi * z[[j]]) * drop (cjj °/ 0 *°/ 0 z [ [ j ] ])) 

> 

diirij <- function (z, h, i, j) { 
cij <- crossprod (h[[i]], h[[j]]) 
cii <- crossprod (h[[i]]) 
czj <- drop (cij %*% z[[j]]) 
czi <- drop (cii z[[i]]) 
rij <- sum (z[[i]] * czj) 
dii <- 3 * rij * outer (czi, czi) 
dii <- dii - outer (czj, czi) 
dii <- dii - outer (czi, czj) 
dii <- dii - rij * cii 
return (dii) 

} 

djjrij <- function (z, h, i, j) { 
cji <- crossprod (hC[j]], h[[i]]) 
cjj <- crossprod (h[[j]]) 
czi <- drop (cji •/,*'/, z [ [i] ]) 
czj <- drop (cjj %*% z[[j]]) 
rij <- sum (z [ [j]] * czi) 
djj <- 3 * rij * outer (czj, czj) 
djj <- djj - outer (czi, czj) 
djj <- djj - outer (czj, czi) 
djj <~ djj - rij * cjj 
return (djj) 

} 

dijrij <- function (z, h, i, j) { 
cij <- crossprod (h[[i]], h[[j]]) 
cii <- crossprod (h[[i]]) 
cjj <- crossprod (h[[j]j) 
szi <- drop (t (cij) %*% z [ [i]]) 
szj <- drop (cij %*'/, z [ [j] ]) 
czj <- drop (cjj '/,*% zCCj]]) 
czi <- drop (cii °/ 0 *°/ 0 z [ [i] ]) 
rij <- sum (z [[i]] * szj) 
dij <- cij + rij * outer (czi, czj) 
dij <- dij - outer (czi, szi) 
dij <- dij - outer (szj, czj) 
return (dij) 

} 
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checkZ <- function (z, h, p, s) { 
fl <- function (z, h, p, s) { 

e <- ff(psplit (z, sapply (h, ncol)), h, p) 
return (e$values[s]) 

> 

fy <- function (z, h, p, s) { 

e <- ff(psplit (z, sapply (h, ncol)), h, p) 

k <- which.max (abs (e$vectors[, s])) 

return (sign (e$vectors[k, s]) * e$vectors[, s]) 

} 

gin <- grad ( 

fl, 

unlist (z), 
side = NULL, 
h = h, 

P = P> 
s = s 

) 

gyn <- jacobian ( 

fy, 

unlist (z), 
side = NULL, 
h = h, 

P = P, 
s = s 

) 

hln <- hessian (fl, 

unlist (z), 
h = h, 

P = P, 
s = s) 

return (list (gin = gin, 
gyn = gyn, 
hln = hln)) 

> 
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